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Abstract 

In this overview paper, we review the basic princi¬ 
ples of the method of space-time conservation ele¬ 
ment and solution element for solving the conser¬ 
vation laws in one and two spatial dimensions. The 
present method is developed on the basis of local and 
global flux conservation in a space-time domain, in 
which space and time are treated in a unified man¬ 
ner. In contrast to the modern upwind schemes, 
the approach here does not use the Riemann solver 
and the reconstruction procedure as the building 
blocks. The drawbacks of the upwind approach, 
such as the difficulty of rationally extending the ID 
scalar approach to systems of equations and par¬ 
ticularly to multiple dimensions is here contrasted 
with the uniformity and ease of generalization of the 
CE/SE ID scalar schemes to systems of equations 
and to multiple spatial dimensions. The assured 
compatibility with the simplest type of unstruc¬ 
tured meshes, and the uniquely simple nonreflecting 
boundary conditions of the present method are also 
discussed. The present approach has yielded high- 
resolution shocks, rarefaction waves, acoustic waves, 
vortices, ZND detonation waves and shock/acoustic 
waves/vortices interactions. Moreover, since no di¬ 
rectional splitting is employed, numerical resolution 
of two-dimensional calculations is comparable to 
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that of the one-dimensional calculations. Some sam¬ 
ple applications displaying the strengths and broad 
applicability of the CE/SE method are reviewed. 


1 Introduction 

This paper provides an introductory overview of a 
new numerical framework for solving conservation 
laws, namely, the Method of Space-Time Conserva¬ 
tion Element and Solution Element, or the CE/SE 
method for short. The method has been developed 
in the past few years by Chang and coworkers [1- 
31]. This method was designed from scratch to be a 
coherent framework, and is not an incremental im¬ 
provement of a previously existing method. All the 
numerical schemes obtained by application of this 
method share the attributes of generality, accuracy 
and simplicity. The method has been applied by 
Chang and coworkers to obtain numerical schemes 
for conservation laws in fluid dynamics. However, 
the method could also be applied to other contin¬ 
uum conservation laws, such as those encountered 
in electromagnetics. 

The CE/SE method is distinguished from other 
methods by its very conceptual basis — flux conser¬ 
vation in space-time. The original reports described 
the present method in a thorough and highly de¬ 
tailed fashion. That is, starting from the basic inte¬ 
gral equations, the algebraic details of the discretiza¬ 
tion and the stability and error analyses were pre¬ 
sented in a systematic way. One of the aims of this 
overview is to provide the reader unfamiliar with the 
CE/SE method with a gentler introduction to the 
method. In this overview, following the treatment of 
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[11], the CE/SE method is described in a more intu¬ 
itive manner by focusing on the geometric nature of 
its unique space-time discretization. This will give 
the reader a complete and simple conceptual descrip¬ 
tion of CE/SE algorithms in one and more spatial 
dimensions. For algebraic details of the schemes and 
numerical analysis of their properties, the reader is 
referred to the original papers. The geometric de¬ 
scription clarifies the conceptual and algorithmic dif¬ 
ferences between the CE/SE method and traditional 
numerical methods for conservation laws. This high¬ 
lighting of these differences is another aim of the 
current paper. The third aim of this paper is to re¬ 
view some sample numerical results obtained with 
CE/SE schemes, in order to give the reader some 
indication of the accuracy and broad applicability 
demonstrated already by the method. A fourth aim 
is to provide the reader with a CE/SE bibliography, 
complete as of July 1998. 

To persuade the reader of the advantages of the 
CE/SE approach, the key features and principal suc¬ 
cesses of the method will now be briefly listed. The 
statements will be supported by the exposition in 
later sections of the paper. Some distinguishing fea¬ 
tures of the CE/SE method are : (i) Although the 
differential form is also utilized, the main empha¬ 
sis of the present method is on solving the inte¬ 
gral form of the conservation equations in the space- 
time domain, (ii) The conservation laws are recog¬ 
nized as balances of space-time fluxes exiting space- 
time regions. The present method was therefore de¬ 
signed to conserve space-time flux locally and glob¬ 
ally. Conventional ‘conservative’ schemes, however, 
calculate spatial fluxes at an instant of time; they do 
not recognize the importance of conservative treat¬ 
ment of the temporal evolution. The space-time 
flux perspective allows for truly conservative treat¬ 
ment of moving boundaries and moving meshes, (iii) 
The discretized equations obtained by the present 
method are a faithful discrete counterpart of the 
conservation laws, automatically reproducing the 
key properties of the latter, such as characteristics- 
related properties. For example, for isentropic flows, 
the present method can be used to construct explicit 
solvers that are non-dissipative (neutrally stable) for 
all Courant numbers <1. Also, these solvers are 
two-way time-marching schemes, i.e., each forward- 
time-marching scheme can be inverted to become a 
backward-time-marching scheme, which is also neu¬ 
trally stable. In other words, the marching variables 
at the (n—l)th time level can be determined in terms 
of those at the nth time level. Thus, these solvers 
reproduce the time-reversibility displayed by isen¬ 


tropic initial-value flow problems arising from pure 
convection phenomena, (iv) A staggered space-time 
mesh is used, i.e., the spatial arrangement of the 
mesh nodes at any time level is spatially staggered 
relative to the arrangement at the previous time 
level. This ensures that flow information at each 
interface separating two conservation elements can 
be evaluated without interpolation (e.g., averaging) 
or extrapolation. In particular, no Riemann solver 
is needed in calculating interfacial fluxes, (v) The 
flow solution structure is not calculated through a 
reconstruction procedure. Instead, the gradients of 
flow variables are treated as independent unknowns, 
and they are not influenced by the flow properties in 
neighboring elements at the same time level. This is 
in full compliance with the flow physics of the initial- 
value problem, (vi) For flows in multiple spatial di¬ 
mensions, no directional splitting is employed. The 
two and three-dimensional spatial meshes employed 
by the present method are built from triangles and 
tetrahedrons, respectively. The CE/SE schemes in 
two and three spatial dimensions are naturally for¬ 
mulated on unstructured triangular and tetrahedral 
meshes which are the most easily generated meshes 
for complicated flow boundary geometries, (vii) The 
explicit time-marching CE/SE algorithms for initial- 
value and initial-boundary-value problems are easy 
to vectorize and parallelize, in order to take ad¬ 
vantage of advanced computer architectures, (viii) 
The flux-based specification of the CE/SE schemes 
give rise in a natural fashion to extremely simple 
yet highly effective non-reflecting boundary condi¬ 
tions. This is in contrast to the complexity of non- 
reflecting boundary conditions necessary for tradi¬ 
tional numerical methods. 

Computer programs based on the CE/SE method 
have been developed for calculating flows in one 
and two spatial dimensions. Numerous results 
were obtained [1-31], including various shock tube 
problems, the ZND detonation waves, the implo¬ 
sion and explosion problem, shocks over a forward- 
facing step, acoustic waves, and shock/acoustic 
wave interactions. The method can clearly resolve 
shock/acoustic wave interactions wherein the differ¬ 
ence of the magnitude between acoustic wave and 
shock could be up to six orders. In two-dimensional 
flows, the reflected shock is as crisp as the lead¬ 
ing shock. From the evidence of these results, the 
CE/SE method has proved to be a promising numer¬ 
ical framework for solving fluid dynamics problems. 

The remainder of the paper is organized as fol¬ 
lows. In Sec. 2, we contrast the present perspective 
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of conservation of space-time fluxes with the tra¬ 
ditional focus on conservation of spatial fluxes. A 
new space-time integral form of conservation laws 
will be described. In Sec. 3, we present the CE/SE 
method for solving flow equations in one spatial di¬ 
mension. In Secs. 4 and 5, the extensions of the 
CE/SE method in two and three spatial dimensions, 
without directional-splitting, are illustrated. In Sec. 
6, numerical examples calculated by the present 
method are shown. Sec. 7 has a summary and con¬ 
cluding remarks. 


2 Numerical Methods from 
the Space-Time Perspective 

2.1 Finite-Volume Conservative Dis¬ 
cretization 

The finite-volume method is the traditional numeri¬ 
cal method that is conceptually closest to the present 
method. Therefore, we begin this section with some 
remarks regarding the basis and limitations of up¬ 
wind finite volume discretizations. Historically, the 
solution of conservation laws was first viewed as 
being reducible to the solution of ordinary or par¬ 
tial differential equations, with specified initial and 
boundary conditions. As the development of nu¬ 
merical methods for differential equations matured, 
it was recognized that flux-conservative numerical 
schemes were essential for the accurate computation 
of flowfields with shocks and other embedded discon¬ 
tinuities. This property of conservative differencing 
was due to its mimicking of integral flux balances, 
which are applicable to discontinuities. With the 
early emphasis on steady state solutions, the conven¬ 
tional finite volume methods for simulating conser¬ 
vation laws were focused on spatial flux balance over 
a fixed spatial domain. Modern upwind techniques 
that stem from Godunov’s scheme, and which fol¬ 
low the projection-upwinding-evolution framework 
described by van Leer [39, 40]. A particularly lu¬ 
cid description of the modern upwind approach is 
found in Huynh [41]. 

Consider conservation of any extensive property 
with continuum distribution per unit volume u, and 
spatial flux vector /. The differential form of the 
conservation or transport law, which applies in the 
absence of discontinuities of the solution and its gra¬ 


dient, is 

du - 

— = -V-/, (2.1) 

where V is the spatial gradient operator. It is this 
differential form that is numerically approximated 
by finite-difference techniques. Integrating Eq. (2.1) 
over a fixed spatial volume V, and using Gauss’ di¬ 
vergence theorem to convert the volume integral on 
the right-hand-side to a surface integral, we obtain 



where S(V) is the boundary of V, and dS = da n 
with dcr and n, respectively, being the area and the 
outward unit normal vector of a surface element on 
S'(V’). Eq. (2.2) is the balance at a particular in¬ 
stant of time between the spatial fluxes exiting a 
finite spatial volume, and the time rate of change 
of a quantity within the volume. The traditional 
finite volume methods concentrate on the conserva¬ 
tive evaluation of the right hand side of Eq. (2.2). 
The left hand side of Eq. (2.2) is often discretized by 
a finite difference method, such as the Runge Kutta 
method. Thus, while the spatial fluxes are treated 
in an integral sense, the flux of the solution in the 
time direction is still treated in a differential man¬ 
ner. This means that the balancing of fluxes occurs 
only at discrete instants of time. The solution in be¬ 
tween these instants will not in general display exact 
balance of the numerical fluxes. This allows the nu¬ 
merical solution to “leak away” as it evolves in time. 
Thus energy and other quantities will generally not 
be globally conserved in an evolving unsteady flow. 

If we integrate Eq. (2.2) with respect to time, we 
obtain 



Note that Eq. (2.3) is more fundamental than Eq. 
(2.1), and applies even to discontinous solutions such 
as shock waves and contact discontinuities. For 
smooth solutions, Eq. (2.1) can be obtained from 
Eq. (2.3) by application of Gauss’ divergence theo¬ 
rem to arbitrary volumes. Eq. (2.3) is indeed an in¬ 
tegral form of the flux balances, and it is the starting 
point for some finite-volume discretizations. How¬ 
ever, even this form has a drawback, namely, that 
the spatial volume V is unchanging in time. This 
fact, coupled with a spatial (rather than space-time) 
perspective on flux conservation, means that in the 
finite-volume context, Eq. (2.3) is generally applied 
to a spatial mesh that is the same at successive time 
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levels. The same kind of mesh is used also when the 
left-hand side of Eq. (2.2) is treated by Runge-Kutta 
integration, or other such method. This awkward 
space-time mesh arrangement evolved from the lim¬ 
ited viewpoint of considering conservation of only 
the spatial fluxes. As will be shown next, it leads 
to the necessity of averaging fluxes at interfaces be¬ 
tween cells, which is in turn the source of much com¬ 
plexity and difficulty. 

As shown in Fig. 2.1(a), due to the fixed spatial 
domain assumed in Eq. (2.3), the shape of the space- 
time Conservation Elements (CEs) in one spatial di¬ 
mension must be rectangular. In addition, because 
of the limited spatial perspective of the mesh in most 
discretizations, these elements are generally stacked 
exactly on top of each other in the time direction, 
i.e., no staggering of these elements in time is used. 
If Eq. (2.2) is used instead, the spatial CEs are just 
the horizontal line segments in Fig. 2.1(a). For equa¬ 
tions in two spatial dimensions, as depicted in Fig. 
2.1(b), a conservation element is a uniform-cross- 
section cylinder in space-time, and again no stagger¬ 
ing in time is employed. This arrangement results in 
vertical interfaces extended in the direction of time 
evolution between adjacent CEs. Across these inter¬ 
faces, flow property information from the previous 
time level travels in both directions. Therefore, in 
calculating the interfacial flux which is needed to 
balance fluxes for the CEs, two values of the interfa¬ 
cial flux are obtained, one from each of the two CEs 
which share the interface. To reconcile these two val¬ 
ues and obtain a unique value of the interfacial flux, 
some form of averaging becomes necessary. A sim¬ 
ple arithmetic averaging of fluxes leads to a central- 
difference approximation of the spatial derivatives. 
In modern high-resolution schemes, to better model 
the flow of information, the averaging is accom¬ 
plished by upwind biasing (or a Riemann solver). 
However, this upwind biasing, termed flux-splitting 
for systems of equations, is based on the Method 
of Characteristics, which is valid only for smooth 
solutions and does not apply to solution discontinu¬ 
ities sucli as shocks. Further, for flows in multiple 
spatial dimensions, directional splitting is used to 
implement one-dimensional characteristic flux split¬ 
ting. This practice causes deterioration of numeri¬ 
cal resolution and difficulties in solving conservation 
laws with source terms, because source terms have 
no preferred direction. 

We mention here that besides the need for the 
averaging of fluxes at interfaces, there are other 
sources of complexity and difficulty in the conven¬ 


tional finite-volume approach. One such source is 
that typically the unknowns in the discretized solu¬ 
tion are just the values of the solution at the mesh 
points. At any mesh point, the numerical approxi¬ 
mations to the flow gradients, which are needed for 
higher accuracy, must be recovered using the dis¬ 
cretized solution at the mesh point and the neighbor¬ 
ing mesh points. This procedure is termed projec¬ 
tion or reconstruction, and is normally accomplished 
by fitting a polynomial to the solution at neighboring 
mesh points. In the vicinity of shocks, however, this 
results in spatial oscillations in the solution. This 
problem is addressed through the use of complex 
flux-limiting strategies, such as Total Variation Di¬ 
minishing and Essentially Non-Oscillatory schemes, 
etc., which use some non-general properties of sim¬ 
ple shock waves. 

Schemes that use a staggered space-time mesh, of 
course, date back to such classical schemes as the 
Leapfrog scheme, which latter scheme in actuality 
computes two decoupled staggered-mesh solutions. 
In general, they are not currently favored by the 
CFD community. However, we point out that there 
are some staggered space-time mesh schemes devel¬ 
oped prior to the current work, which continue to 
be researched. Sanders [47] and Sanders and Weiser 
([48], [49]) developed a staggered mesh scheme. 
However, they did not gain the potential benefit of 
avoiding upwinding, because they used the method 
of characteristics as part of their solution procedure. 
Nessyahu and Tadmor [46] also developed a stag¬ 
gered mesh scheme. However, they did not treat 
the flow gradients as independent unknowns at the 
mesh point, and were thus led to reconstructive pro¬ 
cedures. Both the preceding schemes are dissipa¬ 
tive and therefore irreversible. The present method 
was developed independently in a different frame¬ 
work by Chang and coworkers, without knowledge of 
the schemes of Sanders and Weiser, or of Nessyahu 
and Tadmor. 


2.2 Space-Time Integral Form of the 
Governing Equations 

We now present the space-time integral form of the 
governing equations that is the starting point for the 
CE/SE method. For ease of exposition, we will deal 
with the one-dimensional Euler equations. The one¬ 
dimensional unsteady Euler equations of a perfect 
gas can be expressed as 

LJt -f F r = 0, (2*4) 
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where 


= (/*>), (2.5) 

\PE I 

= ( pv 2 +P ), ( 2 . 6 ) 

\ (pE + p)v / 

with p,v,p, and E being the density, velocity, static 
pressure, and specific total energy, respectively. By 
definition, E = v 2 / 2 + e, where e is the specific in¬ 
ternal energy. The equation of state is p = (7 — l)pe 
with 7 being the specific heat ratio. 


arise when the CE/SE method is applied to mov¬ 
ing boundary or moving mesh problems. Also, the 
space-time perspective inherent in Eq. (2.8) frees 
up the researcher to visualize spatial meshes that 
are not stacked directly one on top of another in the 
time direction, but are instead spatially staggered 
relative to previous time levels. 


2.3 An Ideal Numerical Analogue 



As shown in Fig. 2.2, let x x = x , and x 2 = t 
be the coordinates of a two-dimensional Euclidean 
space £ 2 . Then Eq. (2.4) can be expressed as the 
divergence-free conditions 

V *^m — Oj m — 1,2,3, (2.7) 

where V here denotes the space-time gradient oper¬ 
ator (d/dx,d/dt), and h m = (/ m ,ti m ), m = 1 , 2 , 3 , 
are the space-time mass, momentum, and energy 
current-density vectors, respectively. Equation ( 2 . 7 ) 
is valid everywhere in E 2 for continuous and isen- 
tropic flow solutions. For flows with shock waves, 
we must use the more fundamental form of the con¬ 
servation laws: 

<p h m -ds = 0, m= 1,2,3. (2.8) 

JS(R) V ' 

Here S(R) is the boundary of a space-time region 
R, and ds = dan with da and n, respectively, be¬ 
ing the area and the outward unit normal vector 
of a surface element on S(R). Note that (i) because 
hrn -ds is the space-time flux of h m leaving R through 
ds, Eq. J 2 . 8 ) simply states that the total space-time 
flux of h m leaving R through its boundary vanishes; 
(ii) all mathematical operations can be carried out 
as though E 2 were an ordinary two-dimensional Eu¬ 
clidean space; and (iii) Eq. (2.7) is valid only for 
smooth flows, and it can be derived from Eq. (2.8) 
using Gauss’ divergence theorem. Since it places no 
constraint on the shape of the CEs in the space-time 
domain, Eq. (2.8) is more general than the usual 
mathematical statement of a conservation law, as 
typified by Eq. (2.2) or Eq. (2.3). Eq. (2.8) is the 
space-time integral form that is numerically approxi¬ 
mated in the CE/SE method. For the CEs as defined 
in CE/SE schemes to date, Eq. (2.3) suffices. How¬ 
ever, an alternate derivation of ID CE/SE schemes 
requires the use of Eq. (2.8). The form Eq. (2.8) is 
also essential for the oblique cylinder CEs that will 


A smooth solution to the Euler equations, Eq. (2.4), 
has the following important properties: (i) it does 
not dissipate with time; (ii) its value at any point 
(x, t ) has a finite domain of dependence at an earlier 
time; and (iii) it is completely determined by the 
initial data at a given time. In the light of these 
properties, we remark that (i) a solution to a dis¬ 
sipative numerical scheme will dissipate with time; 
(ii) the value of a solution to an implicit scheme at 
any point ( x,t ) depends on all the initial data and 
all the boundary data up to the time t ; and (iii) a 
scheme involving more than two time levels requires 
the specification of the initial data at more than one 
time level. Therefore, we conclude that an ideal nu¬ 
merical analogue to Eq. (2.4) should be neutrally 
stable , explicit , and involving only two time levels . 

By adding an artificial dissipation term to intro¬ 
duce irreversibility, an ideal solver of Eq. (2.4) can 
be extended to model flows with shocks. We want 
to emphasize that the artificial dissipation in an 
ideal numerical method should occur only in shock 
capturing; without added artificial damping, there 
should be no other source of numerical dissipation. 

Furthermore, in an ideal Navier Stokes solver, the 
above guidelines of modeling the Euler equations 
should be applied to the discretization of the con¬ 
vective terms of the Navier Stokes equations. We 
note that, stripped of any added artificial terms, tra¬ 
ditional numerical schemes are in general not free 
from inherent numerical dissipation. For flows at 
large Reynolds numbers, numerical dissipation may 
overwhelm the physical dissipation and cause a com¬ 
plete distortion of the solution. Because an ideal 
analogue of Eq. (2.4) has no numerical dissipation, 
when it is applied to discretize the convective terms 
of the Navier Stokes equations, the Navier-Stokes 
solver has the property that the numerical dissipa¬ 
tion of its solutions approaches zero as the physical 
dissipation approaches zero. 
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3 The CE/SE Method in One 
Spatial Dimension 

By comparison with modern upwind finite-volume 
methods for hyperbolic equations, the following dis¬ 
tinguishing features of the CE/SE method result in 
a simpler and more consistent numerical flow model: 
(i) A space-time discretization is chosen for the flux 
conservation such that there is no Riemann problem 
to be solved at the cell interface, (ii) To avoid impos¬ 
ing predetermined constraints on the flow solutions 
such as monotonicity and TVD, the flow properties 
and their gradients are treated as unknowns. For 
smooth flows, the unknowns are completely deter¬ 
mined by flux conservation, and the resultant nu¬ 
merical procedure can march forward and backward 
m time, (iii) For flows with shocks, an adjustable 
artificial damping is added to the discretized equa¬ 
tions such that the numerical entropy condition is 
satisfied. 


3.1 Preliminaries 

For simplicity, we develop the CE/SE Euler solver 
on a uniform staggered space-time mesh. In Fig. 3.1, 
we illustrate the nodes, denoted by dots (filled cir¬ 
cles), where the unknowns are located. The space 
and time intervals between neighboring nodal lines 
are respectively denoted by Ax/2 and At/2. Note 
that the spatial interval between nodes at a given 
time level is denoted by Ax. The reason for the no¬ 
tation At/2 for the interval between successive time 
levels is to emphasize that on this staggered space- 
time mesh, it takes two half-time-steps to return to 
the original spatial node locations. There is a Solu¬ 
tion Element (SE) associated with each node (J, n). 
Let the SE(j, n) be the interior of the space-time re¬ 
gion bounded by a dashed line depicted in Fig. 3.2. 
It includes a horizontal line segment, a vertical line 
segment, and their immediate neighborhood. Be¬ 
tween SEs, discontinuities are allowed. For the Eu¬ 
ler equations, Eq. (2.4), which have no source terms, 
the actual size of the neighborhood does not matter. 

Within a SE, the flow property vector U and 
the flux vector F are approximated by their dis¬ 
cretized counterparts U* and F*. Since a second- 
order scheme is desired, piecewise linear functions of 
space and time U* and F* are assumed. For (x,t) 
in SE (j, n), we assume that 

U*0M;i,n) = 


u ” + ( u *)? (* - *j) + (U*)" (t - t n ) (3.1) 

and 

F*(x,t;j, n) = 

( p )“ + ( F «)“ (* " *i) + (F0? (t - t n ). (3.2) 

This is the Taylor polynomial of degree one within 
each element, and the expansion coefficients U?, 
etc., are the column matrices of unknown constants 
to be determined for each SE(j, n). F" is the column 
matrix F (which is a function of TJ) evaluated with 
U = U". 

The expansion coefficients (U/)", (F x )? and (F t )" 
in Eqs. (3.1) and (3.2) will be expressed as functions 
of the independent unknowns U? and (U x )" of the 
present scheme as follows: 3 

(Fx); = a; (u,); , (3.3) 

(u*); = -(f x );, (3.4) 

and 

(fo; = a;(uo;. (3.5) 

Here (i) AJ is the Jacobian matrix A = dF/dU 
(which is a function of U) evaluated with U = U?, 
and (ii) Eqs. (3.3), (3.5) and (3.4) are the numerical 
analogues of the chain-rule relations F x = AU, 

Ft = AU* and the differential equation, Eq. (2.4), 
respectively. ^Furthermore, because the space-time 
flux vectors h m = (/m>u m ), m = 1,2,3, we shall 
assume that for m = 1,2,3, 

—* 

j, n)), (3.6) 

where u m and m = 1,2,3, are the components 
of the column matrices U* and F*, respectively. 

At this juncture, note that hereafter, the compo¬ 
nents of the column matrices U;, (U x )”, (U t )”, F 1 ?, 

(F x )? and (F<)", will be denoted by (u m )" (u mx )j, 

(Umt)j, (fm)j, (fmx)j and m = 1,2,3, re- 

spectively. 


3.2 Discretization of Space-Time 
Conservation Laws 

For smooth flows, the calculations of U; and (U x )" 
are determined by requiring fluxes to be conserved 
over space-time Conservation Elements (CEs). As 
depicted in Figs. 3.3(a) and 3.3(b), two CEs, de¬ 
noted by CE_ (j, n) and CE+(j, n), are associated 
with every mesh point ( j , n). A glance over Figs. 3.1, 
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3.2, and 3.3 reveals that the set of CEi^*, n) over all 
mesh points (j,n) do not overlap among themselves 
and can fill the entire space-time computational do¬ 
main. 


For each (j, n), the following discrete analogues 
to the space-time flux conservation, Eq. (2.8), are 
imposed: 


and 


/ 

JS(CE.(j, n )) 

/ 

JS(CE+(j,n)) 



■ ds = 0, 


• ds = 0, 


m = 1,2,3, 


m = 1,2,3. 


(3.7) 


(3.8) 


According to Figs. 3.2, 3.3a, and 3.3b, we have 
the following observations: (i) The edges CB and 
CD of CE_ (j, n ) l ie in SE(j - 1/2, n - 1/2); (ii) 
The edges AB and AD of CE_(i, n) lie in SE(j, n); 
(iii) The edges ED and EF of CE + (j, n) lie in 
SE(i + 1/2, n - 1/2); (iv) The edges AD and AF of 
CE+(j, n) lie in SE(j,n). As a result, with the aid 
of the numerical counterpart of Eq. (2.8), and Eqs. 
(3.1)-(3.6), we conclude that Eq. (3.7) leads to three 
relations (one for each m) involving the independent 

unknowns U?, (U*)?, U &*)]!$, and Eq. 
(3.8) leads to another three relations, involving U”, 

(U x )”, U^ 1 / 2 , and (Uj)^ 1 ^ 2 . Assuming that the 
unknowns at the mesh points (j - 1/2, n - 1/2) and 
0 + 1/2) n — 1/2) are given, the six components of U/ 
and (U x ) j can be determined by the above six rela¬ 
tions. It will be shown next how JJJ and (U*)? can 
be determined from these relations without needing 
to solve any nonlinear algebraic equations. 


Note that the space-time flux of h* leaving 
CE_ (j, n) th rough AD and that leaving CE +(j, n) 
through AD are evaluated using the same unknowns, 
i.e., U” and (U*)?. Thus, these two space-time 
fluxes are each the negative of the other. As a result, 
a combination of Eq. (3.7) and Eq. (3.8) imply that 

icEox/” " &= °' m=I ' 2 ' 3 ’ (3 ' 9) 

where CE(j, n) (see Fig. 3.3(c)) is the union of 
CE-0» and CE+(j,n). Thus, the CE(j, n) are 
also conservation elements for the scheme. Here, as 
explained above, the fluxes of h* m leaving CE(j, n) 
through CD, CB, ED, and EF can be evaluated in 
terms of U”.^ 2 and (U t )^ 2 2 . Further, it can be 

shown that the flux of h* m leaving CE(j, n) through 
BF is simply (u m )JAx. (This would not be true 


on a nonuniform spatial mesh, but a small modifica¬ 
tion of the method, not described here, allows one to 
similarly avoid solving any nonlinear algebraic equa¬ 
tions.) Thus, Eq. (3.9), which is a combination of 
Eqs. (3.7) and (3.8), implies that U? can be deter¬ 
mined explicitly in terms of U n and (U*)^ 2 . 

After obtaining U?, F” and A” can be determined 
because they are functions of U? only. As a re¬ 
sult, by applying either Eq. (3.7) or Eq. (3.8) (only 
one of these two equations is independent after Eq. 
(3.9) is used), one can obtain a system of three linear 
equations with the three unknowns being the three 
components of (U*)’?. In other words, (U x )j? can be 

determined in terms of Uj?, U””^ 2 , and (U*)”"* 72 
by solving either Eq. (3.7) or Eq. (3.8). 3 

3.3 Special Features of the CE/SE 
Method for Isentropic Flows 

It has been shown by numerical experiments that 
the present scheme is neutrally stable in the inte¬ 
rior of the computational domain up to at least a 
thousand time steps when the Courant number does 
not exceed unity. As a matter of fact, by using an 
analysis similar to that given at the end of Sec. 6 
in [6], one can show that the linearized form of the 
present numerical analogue is neutrally stable when 
the Courant number does not exceed unity. Thus, 
the scheme described above is non-dissipative when 
the Courant number does not exceed unity. This has 
been shown in [2] to result from the fact that the nu¬ 
merical scheme shares with the Euler equations the 
property of being invariant under space-time inver¬ 
sion. The present scheme also meets the require¬ 
ments of an ideal numerical analogue set forth in 
Sec. 2.3, i.e., it should be neutrally stable, explicit, 
and involving only two time levels. 

From the previous development, it can be deduced 
that the numerical scheme observes a global conser¬ 
vation condition that is a direct result of Eq. (3.9), 
i.e., for any space-time region that is the union of 
any combination of the CEs of the type depicted in 
Fig. 3.3(c), the total flux of , m — 1,2,3, leaving 
its boundary vanishes. Thus, the CE/SE scheme en¬ 
forces global conservation of the numerical fluxes in 
both time and space. 

The fluxes are uniquely defined at the boundaries 
of each CE. Each segment of the boundary of a CE 
lies in one and only one SE. This is made possible by 
the spatial staggering of the mesh points at one half- 
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time-level relative to those at the previous half-time- 
level. The conception of this staggered arrangement 
was in turn facilitated by the unified space-time 
viewpoint employed in the CE/SE method. This 
feature of unique definition of the interfacial fluxes 
allows us to avoid the need for upwind or other aver¬ 
aging of the fluxes at the interfaces. Thus the need 
for an exact or approximate Riemann solver does 
not arise. The flow of information between SEs is 
contained within the conservation law itself. 


from Eq. (3.9), the former is incorporated into the 
shock-capturing scheme without modification. As 
a result, given the same and (U^)^ 2 , 

the shock-capturing scheme shares with the non- 
dissipative scheme the same zero-order terms on the 
right sides of Eqs. (3.1) and (3.2). In addition, the 
shock-capturing scheme also observes a global con¬ 
servation condition that is also a direct result of Ea 
(3.9). 4 * 


The CE/SE scheme does not require a flow gra¬ 
dient reconstruction strategy, for the simple reason 
that the flow gradients are themselves treated as 
unknowns, and are determined, together with the 
flow solution itself, from the discretized conservation 
equations, Eqs. (3.8) and (3.9). No unnecessary as¬ 
sumption of smoothness of the solution between SEs 
is made in the Euler isentropic-flow scheme, when 
calculating either the flow variables or the flow gra¬ 
dients. While conventional finite-volume techniques 
also allow for the flow variables to be discontinuous, 
by using conservation to link cells, they do assume 
continuity of the flow when reconstructing the gradi¬ 
ent and when using upwind biasing to calculate the 
interfacial fluxes. Contrarily, in the CE/SE scheme, 
no such assumptions are required, and consequently 
no complex flux-limiter strategy is needed. In some 
finite-element schemes, even more constraining as¬ 
sumptions of solution smoothness are made; the so¬ 
lution itself is taken to be continuous between cells 
by choosing to place single-valued solution nodes on 
the interfaces of cells. 


3.4 The Shock-Capturing Scheme 


The above marching scheme for isentropic flows can 
be expressed as 




H 


(jjn-l/2 yn-1/2 .n-1/2 m xn-1/2^ /« 


and 


(u,)? = 

<?,);;$) <3.u) 


Here H and H x are column-matrix functions. Their 
explicit forms can be obtained from Eqs. (5.20)- 
(5.29) in [6] with the assumption that the viscosity 
t* — 0- In the construction of the shock-capturing 
scheme, the local conservation condition, Eq. (3.9), 
is again assumed. Because Eq. (3.10) follows directly 


The shock-capturing scheme is obtained by mod¬ 
ifying Eq. (3.11). To proceed, let 


— t tW — 1/2 

j± 1/2 


)j=tl/2 — Jo + 


At 


(u.j; 


n —1/2 
± 1/2 ■ 


(3.12) 


i.e., (U ))*-)- 1/2 i s a first-order Taylor’s approximation 
of U at the point (j ± 1/2, n). Thus, 


(U*)7 = 


(UT 


j + 1/2 


(u')”_ 


1/2 


Ax 


(3.13) 


is a central-difference approximation of U x at the 
mesh point ( j,n ). In the shock-capturing scheme, 
Eq. (3.11) is replaced by 


( u *)i = (1 - 2e)(H,); + 2e(U')”, (3.14) 

where (H*)? denotes the expression on the right 
side of Eq. (3.11), and e is a real number. Note 
that (U£)” is defined in terms of a central-difference 
approximation. Generally, numerical dissipation is 
introduced as a result of using such an approxi¬ 
mation. On the other hand, (H r )? represents the 
solution from a non-dissipative scheme. The right 
side of Eq. (3.14) is a weighted averaged of (Ha,) 1 ? 
and (U£)? with the weight factor of 1 - 2c and 2c, 
respectively. Therefore, one may heuristically con¬ 
clude that the numerical dissipation associated with 
the shock-capturing scheme can be increased by in¬ 
creasing the value of c. This conclusion is verified by 
numerical experiments. As shown in [6], the stabil¬ 
ity domain of the shock-capturing scheme is defined 
by 


CFL < 1 and 0 < e < 1, 

(3.15) 

where CFL is the maximum Courant number. Note 
that Eq. (3.14) can also be expressed as 

(Ux)” = (Hx)" + 2c(DU)", 

(3.16) 

or 


(Ux)” = (U‘) + (2c - i)(mj)?, 

(3.17) 

where 


(z?u); = (u«);-(h,);. 

(3.18) 
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According to Eq. (3.16), (U X )J for the shock¬ 
capturing scheme is the sum of the non-dissipative 
term (H*)? and the dissipative term 2e(DV)J. The 
latter provides the necessary entropy-increasing con¬ 
dition within the stability domain defined by Eq. 
(3.15). Also it is seen from Eqs. (3.16) and (3.17) 
that (U,)J reduces to (H x )f and (U*)? in the cases 
of c = 0 and e = 0.5, respectively. 


The shock-capturing scheme described above gen¬ 
erally can capture shocks with high resolution and 
without generating substantial numerical oscilla¬ 
tions near shock if 0.3 < e < 0.8. To further damp 
out these oscillations, (U£)?in Eq. (3.17) (which is 
equivalent to Eq. (3.14)) can be modified using a 
simple slope-limiting procedure [ 6 ]. Let 




(3.19) 


Because (U ')? ±1 ^ 2 and UJ are the numerical ana¬ 
logues of U at the mesh points (j± 1/2, n) and (jf, n), 
respectively, (U*+)*? and (U*_)" are two numerical 
analogues of U* at the mesh point (j, n), with one 
being evaluated from the right and another from the 
left. It follows from Eqs. (3.13) and (3.19) that 


( U *); - 5 [(U,+)7 + (U,_)?] , (3.20) 


( u £)? is the simple average of ( 11 *+)*? and 
(U.-)?. A nonlinear weighting function is defined 
as 


W(x~, a?+; a) 


\x+\ a x- 4 - |ar^| Qf x+ 
jaj+l" + x- a 


(3.21) 


where ar+, and a are real variables with 

l*+l + I*-1 > 0 and ot > 0. Note that (i) 
W(x -, # +; 0 ) is the simple average of and a?+, and 
(ii) W(x~, a?+; 1 ) and W(x-,x+ \ 2 ) are used in the 
slope-limiters proposed by van Leer [40] and van Al- 
bada et al. [50], respectively. 


Recall that (u mx ±y- denotes the m-th component 
of (U*±)?. Let 

«x)j = W [(u mr _ )?, (u m *+)", a] (3.22) 

Then, as shown in [ 6 ], numerical oscillations near 
shocks can be suppressed very efficiently if (U£)” in 
Eq. (3.17) is replaced by (UJ?)?, i.e., the column 
matrix formed by m = 1 , 2 , 3 , if a > 1 . 


3.5 Concluding Remarks 


have been developed for the scalar advection equa¬ 
tion (with and without controllable added numerical 
dissipation), and for the scalar advection-diffusion 
equation. Corresponding schemes have been formu¬ 
lated for the Euler and Navier-Stokes equations. The 
Euler schemes were described in this section. Im¬ 
plicit schemes have also been formulated for the ID 
advection-diffusion equation. The various schemes 
for the scalar equations have been thoroughly anal¬ 
ysed for their accuracy and stability properties, see 
[ 6 ] and [ 8 ]. Surprisingly, various limiting cases of the 
ID scalar schemes turn out to have amplification 
tictors matching those of several classical numeri¬ 
cal schemes, viz., the Leapfrog, Lax, DuFort-Frankel 
and Crank-Nicolson schemes. 

The numerical boundary conditions have not been 
discussed in this paper. The boundary conditions 
used to date have been both simple and effective. 
The flux-based nature of the method allows, in par¬ 
ticular, the use of extremely simple yet robust non- 
reflecting boundary conditions [12]. The efficacy of 
these non-reflecting boundary conditions is demon¬ 
strated in some of the numerical examples in this 
paper. Further research is under way on boundary 
conditions for the CE/SE schemes. 


4 The 2D Euler Solvers 


In Sec. 3, it was established that in the ID case 
there were only two sets of independent marching 
variables, i.e., (i) (w m )^ , tti = 1 , 2 , 3 , and (ii) (w m;r ) •, 
m = 1, 2, 3, at each mesh point (j, n), if Eqs. (3.1)- 
(3.6) are assumed. As a result, it requires two sets of 
conservation conditions, i.e., Eqs. (3.7) and (3.8) to 


construct the ID non-dissipative Euler scheme. As 
a prerequisite to Eqs. (3.7) and (3.8), two CEs, i.e., 
CE-C;, n) and CE+(y, n) are defined for each mesh 
point (j,n). 


The 2D CE/SE non-dissipative Euler solver [5, 7 ] 
was constructed using the same set of design princi¬ 
ples that was used to construct its ID counterpart. 
The differences between them stem entirely from the 
fact that there is one more spatial dimension to be 
considered in the 2D solver. In this section, only the 
basic geometric structures of the 2D solver will be 
described. For other details, the reader is referred 
to [5]. 


We conclude this section on the CE/SE method in 
ID by mentioning that ID CE/SE explicit solvers 
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The 2D unsteady Euler equations of a perfect gas 
[5, 7] consist of four independent equations, m = 1 , 



2, 3, 4, instead of the three equations applicable to 
ID flow. Also, in the 2D case, there are two spatial 
components of the gradient of each u m (i.e., u mx 
and Ufnyj where x and y are Cartesian coordinates 
for the 2D space). This is in contrast to the ID 
case, in which, for each u m , there is only one spatial 
component of the gradient of u m (i.e., u mx ). 

In the development of the 2D non-dissipative Eu¬ 
ler solver and its extensions [5], a set of equations 
that is a natural 2D extension of Eqs. (3.1)—(3.6) is 
assumed. As a result, there are three sets of indepen¬ 
dent marching variables at each mesh point (j, l k, n) 
(see Figs. 4.1 and 4.2 for the locations of the mesh 
points. The reader is referred to [5, 7] for the defini¬ 
tions of the spatial mesh indices j and k of the uni¬ 
form structured triangular mesh of Fig. 4.1). They 

are ( Um )j,k> ( u mx)j } k and (u my )^ k , m = 1, 2, 3, 4. 
It follows that it requires three sets of conservation 
conditions (each set comprises four conditions, cor¬ 
responding to m = 1, 2, 3, 4) at each mesh point to 
construct the 2D non-dissipative solver. Therefore, 
as a prerequisite, one must define three conservation 
elements for each mesh point. The construction of 
these CEs, which is the most intriguing part of the 
development of the 2D CE/SE Euler solver, will be 
described in what immediately follows. 

Consider a spatial domain formed by congruent 
triangles (see Fig. 4.1). The center of each triangle 
is marked by either an empty circle or a filled circle. 
The distribution of these empty and filled circles is 
such that if the center of a triangle is marked by a 
filled (empty) circle, then the centers of the three 
neighboring triangles with which the first triangle 
shares a side are marked by empty (filled) circles. 
As an example, point G, the center of the triangle 
A BDF, is marked by a filled circle while points A, C 
and E, the centers of the triangles A BFM, ABJD 
and A DLF , respectively, are marked by empty cir¬ 
cles. These centers are the spatial projections of the 
space-time mesh points used in the 2D solver [5, 7]. 

To specify the exact locations of the mesh points 
in space-time, one must also specify their temporal 
coordinates. In the 2D CE/SE development, again 
we assume that the mesh points are located at the 
time levels n = 0, ±1/2, ±1, ±3/2, ..., with t = 
n A t at the nth time level. Furthermore, we assume 
that the spatial projections of the mesh points at a 
whole-integer (half-integer) time level are the points 
marked by empty (filled) circles in Fig. 4.1. 

Let the triangles depicted in Fig. 4.1 lie on the 
time level n = 0. Then those points marked by 


empty circles are the mesh points at this time level. 
On the other hand, those points marked by filled 
circles are not the mesh points at the time level n — 
0. They are only the spatial projections of the mesh 
points at half-integer time levels. Thus the 2D space- 
time mesh is, as in the ID case, a staggered mesh. 

Points A, C and E , which are depicted in Figs. 4.1 
and 4.2, are three mesh points at the time level n = 
0. Point G', which is depicted in Fig. 4.2, is a mesh 
point at the time level n = 1/2. Its spatial projection 
at the time level n = 0 is point G. Because point G 
is not a mesh point, it is not marked by a filled circle 
in the space-time plot given in Fig. 4.2. Hereafter, 
only a mesh point, e.g., point G 7 , will be marked by 
a filled or empty circle in a space-time plot. 

The conservation elements associated with point 
G f are defined to be the space-time quadrilateral 
cylinders GFABG'F'A'B ', GBCDG 1 B'C'D' , and 
GDEFG'D'E'F f that are depicted in Fig. 4.2. Here 
(i) points B , D and F are the vertices of the trian¬ 
gle with point G being its center (centroid) (see also 
Fig. 4.1), and (ii) points A', 5', C 7 , D 7 , E f and F 1 
are on the time level n = 1/2 with their spatial pro¬ 
jections on the time level n = 0 being points A, B , 
G, D, E and F , respectively. 

Recall that, in the development of the ID non- 
dissipative Euler solver, a pair of diagonally oppo¬ 
site vertices of each CE ± (j,n) (see Figs. 3.3(a) and 
(b)) are assigned as mesh points. Furthermore, the 
boundary of each CE±(j, n) is a subset of the union 
of the SEs associated with the two diagonally oppo¬ 
site mesh points of this CE. In the 2D development, 
as seen from Figs. 4.2, two diagonally opposite ver¬ 
tices of each CE are also assigned as mesh points. 
In the following, we shall define the SEs such that 
even in the 2D case, the boundary of a CE is again 
a subset of the union of the SEs associated with the 
two diagonally opposite mesh points of this CE. 

As an example, the SE associated with 
point G 7 is depicted in Fig. 4.3. It is the 
union of three vertical rectangles (i.e., G ,f B n BG, 
G n D n DG and G N F /, FG) i a horizontal hexagon (i.e., 
A!B ! C*D f E‘F') and their immediate neighborhood. 
Note that points G 77 , B n , D n and F n are on the time 
level n = 1 and their spatial projections on the time 
level n = 0 are points G, B, D and 1 F, respectively. 
The definition of the SE of any mesh point is similar 
to the definition of the SE of the point G 7 . Note that 
on the uniform structured mesh shown, the SEs and 
CEs at the whole integer time levels can be seen to 
be congruent to the SEs and CEs respectively at the 
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half integer time levels, by a rotation of 180 degrees 
about the time axis followed by spatial and temporal 
translations. 

As depicted in Fig. 4.2, one of the CEs associated 
with point G 7 is the space-time quadrilateral cylin¬ 
der GFABG'F'A* B f . Among the vertices of this 
CE, only points A and G 7 are mesh points. From 
Figs. 4.3, it is seen that (i) three of the faces of this 
CE, i.e., GBG'B', GG*F f F and G , F f A , B / are sub¬ 
sets of the SE of point G 7 , and (ii) the other three 
faces, i.e., AA'F'F , ABB 1 A! and ABGF are sub¬ 
sets of the SE of point j4. As a result, by assuming 
that the flux of each h* m (m = 1, 2, 3, 4) is con¬ 
served over this CE, one can impose four conditions 
involving only the independent marching variables 
at the mesh points A and G 7 . Similarly, by using 
the flux conservation conditions over the other two 
CEs associated with point G 7 , one can obtain eight 
other conditions that relate the independent march¬ 
ing variables at the mesh points G 7 , G and E. Using 
the above 12 conditions, the 12 independent march¬ 
ing variables, i.e., u m , u mx and m — 1, 2, 3, 
4, at the mesh point G 7 can be determined in terms 
of the independent marching variables at the mesh 
points A, C and E. By considering the mesh point 
G 7 as a typical mesh point, the reader can under¬ 
stand how the 2D non-dissipative Euler solver was 
constructed [5, 7]. 

The non-dissipative Euler solver is only one of sev¬ 
eral 2D solvers described in [5, 7]. The latter docu¬ 
ment includes the 2D extensions of all but one of the 
ID solvers described in [6]. The only exception is the 
2D extension of the ID Navier-Stokes solver. This is 
under development, and will be dealt with in a sep¬ 
arate paper. Also, because of the similarity in their 
design, each of the 2D extensions shares with its ID 
version virtually the same fundamental characteris¬ 
tics. As an example, the 2D non-dissipative Euler 
solver is neutrally stable, explicit, and involves only 
two time levels during a single time step. It also pre¬ 
serves the forward-backward marching nature and 
the space-time inversion invariance property of the 
2D unsteady Euler equations. These are the same 
properties that characterize the ID non-dissipative 
Euler solver. 

The discussion of the 2D Euler solvers is concluded 
with the following remarks: 

1. Because (i) the spatial geometric structure em¬ 
bedded in the CE/SE 2D space-time mesh is 
constructed from triangles, and (ii) triangles 


are the simplest polygon in the 2D space, the 
CE/SE solvers described in [5] can easily be 
modified and extended to solve flow problems 
with complex geometries using unstructured tri¬ 
angular spatial meshes[31]. 

2. Several 2D CE/SE solvers using nonuniform 
mesh have been developed [22-29]. Some of the 
numerical results generated with these solvers 
will be presented in Sec. 6. 

3. The extension of the non-dissipative Euler 
solver to become a shock-capturing scheme by 
the addition of controllable numerical dissipa¬ 
tion is a straightforward extension of the devel¬ 
opment in Eqs. (3.10) - (3.22). 

4. The non-reflecting boundary conditions men¬ 
tioned in the previous section have also easily 
been extended (e.g., [13, 30]) to the 2D case, 
and have proven to be just as effective as in ID 

5 The Basis for a 3D Euler 
Solver 

We indicate the discretization of space-time which 
forms the basis for a 3D Euler solver currently under 
development. The extension of the CE/SE method 
to three spatial dimensions follows reasoning similar 
to that used when extending the ID solver to the 
2D case (see the previous section). In the 3D case, 
the unsteady Euler equations of a perfect gas con¬ 
sist of five independent equations, m = l,2,3,4,5. 
There are three spatial components of the gradient 
of each u m (i.e., u mx , u my and u mz , where x , y and z 
are Cartesian coordinates for the 3D space). When 
piecewise linear variation with space and time are as¬ 
sumed for the numerical solution, as is done in the 
ID and 2D cases, and after the differential equation 
is assumed valid at each mesh point, there remain 
four sets of independent marching variables at each 
mesh point. It follows that four sets of conserva¬ 
tion conditions are required at each mesh point to 
construct the non-dissipative 3D solver. Hence, four 
conservation elements must be defined for each mesh 
point. Just as a triangle was the polygon sharing its 
bounding edges with three neighbors, so a tetrahe¬ 
dron is the polyhedron sharing its bounding surfaces 
with four neighbors. 

In the 2D case, referring to Figs. 4.1 and 4.2, 
GFAB , GBCD and GDEF are the spatial projec¬ 
tions of the CEs associated with G 7 . The CEs in 
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the 3D case can be constructed in analogous fash¬ 
ion. Consider the tetrahedron ABCD with centroid 
G, and the tetrahedron ABCP with centroid H , 
depicted in Fig. 5.1. They share the face ABC. 
The polyhedron GABCH is then defined as the spa¬ 
tial projection of a CE associated with a point G / . 
The CE is thus a right cylinder in space-time, with 
GABCH as its spatial base. The point G is the spa¬ 
tial image of the mesh point G\ which is displaced 
temporally from G by half a time step. 

In similar fashion, three additional CEs associ¬ 
ated with the mesh point G' can be constructed 
by considering in turn three tetrahedra that share 
with ABCD one of its other three faces. Thus the 
numerical solution at G f can be determined from a 
knowledge of the solution at the four mesh points 
(one of which is H) which are the centroids of the 
tetrahedra sharing a face with ABCD. This forms 
the basis of a non-dissipative 3D Euler solver. 

Just as the structured mesh of Fig. 4.1 is obtain¬ 
able by sectioning the parallelograms of Fig. 4.1 into 
triangles, so it is possible to construct a structured 
mesh of tetrahedra by sectioning a mesh of paral¬ 
lelepipeds. Details of the construction will be given 
in a future paper. Again, the extension to a space- 
time mesh built from an unstructured tetrahedral 
spatial mesh is simple. 


6 Computational Examples 

6.1 Shock Tube with Non-Reflecting 
Computational Boundaries 

The CE/SE computational results are presented for 
an extended Sod’s shock tube problem [51], in which 
the shock tube problem is extended by imposing a 
non-reflecting boundary condition at each end of the 
computational domain. The challenge of the non- 
reflective boundary condition is no less difficult than 
that of capturing the shock and the contact disconti¬ 
nuity. First, the flow under consideration is subsonic 
throughout and the treatment of the non-reflecting 
boundary condition for a subsonic flow is more dif¬ 
ficult than that for a supersonic flow. Second, this 
difficulty is exacerbated by the existence of a shock 
and a contact discontinuity, which must be allowed 
to exit the domain without reflection. 

Flow of an ideal gas with specific-heat ratio 7 = 


1.4 is considered in an infinite shock-tube. The ini¬ 
tial condition, at time t = 0, is ( p, = (l 0 1) 
if * < 0, and (p,v,p) = (0.125,0,0.1) if x > 0. 
Here, p, v,p denote the density, velocity, and pres¬ 
sure of the fluid, respectively. A uniform space-time 
mesh with Ax = 0.01 and At = 0.004 (correspond¬ 
ing to a maximum Courant number of about 0.88) 
is used over the computational domain defined by 
—0.5 < x < 0.5 and t > 0. The settings e = 0.5 and 
cv = 1 are used for the artificial dissipation parame¬ 
ters. Note that the results are obtained without the 
need of any local mesh-refinement techniques or any 
time-step tuning. 

The non-reflecting boundary conditions used are 
(i) u ? = U and (U )” = if 

is a mesh point on the right boundary, and (ii) 

u ” = (u*>7 = (».);;$ » o,n)» 

a mesh point on the left boundary. The reasons 
why such trivial extrapolations can serve so well 
as non-reflecting boundary conditions in the CE/SE 
method are explained in a separate paper [12]. 

Figs. 6.1-6.3 show the numerical solution (tri¬ 
angular data points) compared with the analyti¬ 
cal solution (unbroken line) at three different times, 
namely, t = 0.2, 0.4 and 0.6. It is seen that ex¬ 
cellent agreement is obtained between the numerical 
results and the analytical solution. In particular, 
as seen in Fig. 6.1, the shock wave discontinuity is 
resolved almost within one mesh interval and the 
contact discontinuity is resolved in four mesh inter¬ 
vals. Fig. 6.2 shows that by t = 0.4, the numerically 
computed shock wave has passed cleanly through the 
right boundary, with no spurious reflections. Simi¬ 
larly, Fig. 6.3 shows that by t = 0.6, the contact 
discontinuity has passed through the right bound¬ 
ary, while the expansion region has partially passed 
through the left boundary. Agreement with the ex¬ 
act solution continues to be excellent. 


6.2 Convection-Diffusion Examples 

The CE/SE computations described in this sub¬ 
section were originally presented in [ 8 ], where an 
implicit CE/SE solver for the convection-diffusion 
equation u t + o.u x — pu xx = 0 (/i > 0 ) was devel¬ 
oped. The solver, termed the a-/i(/l) solver, is an 
extension of the a scheme, which is the CE/SE solver 
for u t -f au x = 0 . The examples below help show 
that the scheme is accurate over the whole Reynolds 
number range, from pure diffusion to convection- 
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dominated solutions. 

Pure Diffusion. We consider a special case of 
the convection-diffusion equation with a = 0 and 
H = 1, in the domain 0 < x < 1 and t > 0. The 
initial/boundary conditions completing the prob¬ 
lem specification are (i) u( 0 ,f) = u(l,t) = 0 for 
t > 0 , (ii) ti(*, 0 ) = 2x for 0 < x < 0 . 5 , and (iii) 
u(a?,0) = 2(1 — x) for 0.5 < x < 1 . The solution 
u(x,t) exhibits the diffusive decay of the initial saw¬ 
tooth shape. An exact series solution is available, see 
for e.g. p.15 of [52]. For the CE/SE computation, 
uniform mesh intervals Ax = 0.02 and At = 0.005 
are used. Fig. 6.4 shows the time-slice at t = 0.05, 
comparing numerical and exact solutions, and also 
showing the error scaled with the peak exact value 
at that time level. The maximum error magnitude 
is seen to be about 0.5% of the peak solution value. 
At t = 1 (not shown), when the peak solution value 
has dwindled to about 4 x 10“ 5 , the maximum er¬ 
ror magnitude is about 0.15% of the peak solution 
value. 

Boundary Layer, Re = 100. We next consider 
the problem defined for the convection-diffusion 
equation in the domain 0 < x < 1 and t > 0 by the 
conditions (i) u( 0 ,f) = 0 for t > 0 , (ii) u(l,f) = 1 
for t > 0, and (iii) u(ar,0) = x for 0 < x < 1. The 
‘steady-state* or time-asymptotic limit of the solu¬ 
tion is u(x, oo) = [exp (ax/fi) - 1 ] / [exp(a/^) - 1 ]. 
The case a = 1 , /i = 0.01 (i.e. Re = 100) is consid¬ 
ered, which leads to a steady-state boundary layer 
at x = 1 . Uniform mesh intervals Ax = 0.0025 and 
At = 0.002 are used, so that the Courant number is 
0.8. Fig. 6.5 shows the computed and exact steady- 
state limits, together with the error. The boundary 
layer is seen to be well resolved, with the maximum 
magnitude of the error being about 1 % of the solu¬ 
tion peak. 


6.3 Unstable ZND Detonation Wave 

1 The Piston Problems 

We consider ID and 2D ZND detonation wave 
problems. The governing equations consist of the 
Euler equations together with an equation for the 
reactant concentration, with a stiff source term in 
it. For details of the treatment of source terms in 
the present CE/SE method, we refer the reader to 
[18]. 

The initial condition of the present calculation is 


a long tube filled with reactant with a piston on one 
end moving at a constant speed into the quiescent 
reactant. Here, we use the piston face as the origin 
of the coordinate system. In this coordinate frame, 
reactant is charged into a closed-end tube at a con¬ 
stant speed. A shock wave is reflected on the closed 
end to ignite the reactant. 

The parameters of the flow field in the present 
calculation are set as chemical potential g 0 = 50, ac¬ 
tivation energy E + = 50, specific heat ratio 7 = 1.2, 
and the over drive coefficient / = 1.6. According 
to the classical theory for detonation instability, the 
detonation wave becomes unstable and a longitudi¬ 
nal wave bouncing between the piston and the shock 
front should be observed. In this calculation, only 5 
mesh nodes are used in the each half-reaction zone. 
In Fig. 6 . 6 (a), we show the temporal evolution of the 
pressure level at the shock front. The first pressure 
jump in the figure is caused by the start-up process 
of the pushing piston. After the first pressure jump, 
the flow field settles down and the instability waves 
gradually built up. After 40 time units, a remark¬ 
able instability wave occurs. In about 60 time unit, 
there are about 8 pressure peaks. This numerical 
solution is in excellent agreement with the results 
reported by Fickett and Wood [53]. 

2 Numerical structure for two-dimensional 
detonations 

Unstable detonation waves obtained with the 2D 
Euler equations plus a ZND equation are computed 
by S.T. Yu and S.J. Park, in work yet to be pub¬ 
lished. Fig. 6 . 6 (b) is a Schlieren-type image of deto¬ 
nation waves, especially plotted two periods for pres¬ 
sure. A detonation is traveling from the top to the 
bottom and the flowfield is composed of: (i) the qui¬ 
escent state of the reactant before the shock, (ii) a 
von Neumann spike with finite rate reaction, and 
(iii) the equilibrium state after the reaction zone. 
The two-dimensional cellular structure of detona¬ 
tion waves is of concern: that is, the wave patterns 
that arise when the flow parameters are chosen such 
that the flow field is unstable. The parameters of 
the flow field in the present calculation are set as 
qo = 50, E+ = 50, 7 = 1.2, and the over drive coef¬ 
ficient / = 1 . 6 . According to the classical theory for 
detonation instability, an unstable detonation wave 
should be obtained with these parameters. 

In this calculation, 20 mesh nodes are used in 
the half reaction zone (260 mesh cells covered the 
width of the channel). The calculation of two- 
dimensional detonation wave is initiated by placing a 
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one-dimensional solution on a two-dimensional mesh 
with periodic boundary conditions along the left and 
right of the flow. Without artificial perturbation, 
the shock front quickly becomes unstable under flow 
conditions receptive for instability. Crisp cellular 
structure is observed after transverse waves cause 
the shock front to be curved. This figure shows that 
regular mach stem cell structure exists between two 
curved strong incident shocks and that peak pressure 
zone is around triple points in the cellular structure 
for unstable detonations. 

The results compared favorably with the previ¬ 
ously reported data. These results show that the 
CE/SE method is a very accurate method for direct 
calculations of propagating detonation waves. Fig¬ 
ure 6.6(b) is obtained at t= 12.5sec. 


6.4 Diffraction of a Shock Wave 
around a Wedge 

According to the experimental shadowgraph results 
shown in [54], when a plane shock wave of M s =1.3 
is moving over the beginning of a finite wedge of 
semi-vertex angle 6 = 26.565°, an ordinary Mach 
reflection is generated. As the shock wave passes 
the base, the flow separates to form vortex sheets at 
the sharp corners. Further interaction produces an 
increasingly elaborate pattern of shock waves, slip 
lines and vortices. 

As reported in [22, 26], this flowfield is simulated 
using the CE/SE Euler solver. By virtue of the 
symmetry in the solution, attention is restricted to 
the upper half of the domain. The extent of the 
computational domain is set based on an estima¬ 
tion from Fig. 522 in [54]. The shock wave is at 
x = —0.5 at t = 0. The numerical boundary con¬ 
dition imposed on the vertical wall of the wedge is 
described in [25]. Numerical solutions at eight time 
levels (t = 0.725, 0.9075, 1.2125, 1.55, 1.825, 2.1375, 
2.4875, and 2.9475), obtained by using two subdo¬ 
mains with 321x89 and 209x34 mesh points, and 
with At = 0.0025, are shown in Fig. 6.7. It should 
be pointed out that the upper and lower walls of the 
channel shown in the shadowgraphs of [54] are actu¬ 
ally further apart than the top and bottom edges of 
the shadowgraphs. Therefore some flow phenomena 
that are seen in Fig. 6.7, in the region near the upper 
wall, are lost in shadowgraphs, especially at the 4th 
time level. Comparisons of the computed solutions 
with experimental pictures of [54] have shown an ex¬ 
cellent agreement in general flow features except for 


those phenomena induced by the effect of viscosity. 

The shock waves, slip lines and vortices are captured 
very well. 


6.5 Implosion/Explosion of Polygo¬ 
nal Shock Waves in a Box 

The 2-D CE/SE Euler solver has been used in [27] 
to solve a problem studied in [55], concerning the 
implosion/explosion of a polygonal shock wave in a 
square box. In addition to the early stage of the 
implosion/explosion process, the later development 
of the process, which was not studied in [55], is also 
simulated in [27], The computation further demon¬ 
strates the robustness of the CE/SE Euler scheme 
in handling discontinuous flows. 

A uniformly distributed 241x241 grid is utilized 
in the computational domain, which is a square de¬ 
fined by —2 < a? < 2 and -2 < y < 2. The initial 
shock wave configuration is a polygon, the geomet¬ 
ric center of which coincides with that of the square. 
Inside the polygon is a low pressure region, with a 
pressure ratio of 10 across the shock. The radius of 
the circumscribed circle of the polygon is selected to 
be 0.8\/^ for all shapes of the polygon. In the nu¬ 
merical scheme, the two parameters c and & are set 
to be 0.5 and 1 respectively, everywhere in the com¬ 
putational domain for all cases, and the maximum 
Courant number is always kept at a value of 0.9. 

In one set of computations, the early flowfield is 
studied for polygonal shock waves with initial shapes 
of an equilateral triangle, a square, and a pentagon. 
The density contour plots at different time levels are 
shown in Fig. 6.8. Wave patterns similar to those 
captured in Figs. 1-5 of [55] using a TVD method 
on a 359x359 grid are clearly shown in the CE/SE 
solutions, displaying detailed features such as Mach 
stems and the newly-developed smaller polygons. 

In another computation, the implosion/explosion 
of a hexagonal shock wave is simulated until the re¬ 
implosion of the shock wave is observed in the box. 
More complex flow phenomena can be seen in the 
density contour plots of Fig. 6.9, including the reflec¬ 
tions of shock waves, shock-shock interaction, and 
shock-contact surface interaction. It is interesting 
to note that the shape of the contact surface cen¬ 
tered at the origin of the box remains unchanged 
even after the passage of shock waves. 
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6.6 Examples from Computational 
Aeroacoustics 


consequent reflection amounts to about 1% of the 
strength of the incident waves. 


The CE/SE computational examples we describe in 
the next two subsections were reported in [16] and 
[15]. The investigations in [16] and [15] found the 
CE/SE Euler scheme to be capable of handling the 
complete spectrum of flows, from small-amplitude 
linear acoustic waves, all the way to nonlinear or 
even discontinuous waves (shocks). Through nu¬ 
merical experiments in computational aeroacoustics, 
the following salient properties of the CE/SE Eu¬ 
ler scheme emerge: (i)The CE/SE scheme possesses 
very low dispersion error and yields high resolution 
results comparable to that of a high order compact 
difference scheme, although nominally the CE/SE 
scheme is only of 2nd order accuracy, (ii) In gen¬ 
eral, the numerical non-reflecting boundary condi¬ 
tion applicable to the CE/SE scheme is genuinely 
multi-dimensional, and can be implemented in a sim¬ 
ple and elegant way without resorting to the com¬ 
plexities of characteristic forms or buffer zones, (iii) 
The CE/SE scheme is both a CFD ( Computa¬ 
tional Fluid Dynamics) and a CAA (Computational 
Aeroacoustics) scheme, capable of handling contin¬ 
uous and discontinuous flows. It thus represents a 
unique numerical technique for flows where sound 
waves and shocks and their interactions are impor¬ 
tant. 

It is well-known that in CAA, the non-reflecting 
boundary condition plays a dominant role in the fi¬ 
nal numerical results. In general, there are three 
ways to impose the non-reflecting boundary condi¬ 
tions, namely, 

(i) to apply 1-D characteristic variables (Riemann 
invariants) in the direction normal to the boundary, 

(ii) to minimize spurious numerical reflections 
from the boundaries by inserting a buffer zone with 
increased numerical damping, 

(iii) to apply an asymptotic analytical solution at 
the boundaries. 

In the new CE/SE scheme, none of the above 
complex treatments of non-reflecting boundary con¬ 
ditions is needed. Instead, one of several possible 
simple non-reflecting boundary conditions is: 

( u mx)j i k =: ( u my)^ t k = 0, TYl = 1,2, 3, 4, 

while (^m)”jfc is defined by zeroth order extrapola¬ 
tion from the interior neighbors. In general, the 


6.7 Multiple interaction of strong 
vortex and shocks 

In the first example, the multiple interaction of a 
strong vortex and shocks is considered. It is demon¬ 
strated as a typical example to show that the CE/SE 
scheme is capable of handling not only .the linear 
waves but also the highly nonlinear waves as well. It 
is well-known that interaction between a vortex and 
a shock may generates acoustic waves. To our knowl¬ 
edge, researches including both experimental and 
numerical computation, are mostly concentrated in 
a single interaction between a vortex and a normal 
shock. In this problem, we intend to show the capa¬ 
bility of the CE/SE method in handling complicated 
wave interactions. Although there are no existing 
experimental results for comparison, our numerical 
result appears consistent with the physical phenom¬ 
ena. A grid of 401 x 101 is employed in this problem 
with Ax = Ay = 1. The inflow boundary condition 
is given as a supersonic flow of Mach 2.9: 

u 0 = 2.9, vo = 0, po = 1/1.4, po = 1 

boundary condition at the top is an inclined flow: 

u* = 2.6193, v t = -0.50632, p* = 1.5282, p t = 1.7000 

The outflow boundary condition is of the non¬ 
reflecting variety and the bottom boundary is a solid 
reflecting wall. Then, a steady oblique shock is 
formed with 29° inclination and reflected at the bot¬ 
tom wall. The flow with shocks is precalculated un¬ 
til a steady state is reached. It is then used as the 
background mean flow for further computation. At 
t = 0, a strong Lamb’s vortex is placed at x = 22, 
2/ = 60. The following is a brief description for a 
stationary Gaussian type of Lamb vortex. In polar 
coordinates, the azimuthal and radial velocities uo 
are given as 


u& 



u r 


= 0 . 


With the prescribed u$ and u r , using the momentum 
and energy equations: 


dp ___pu | 
dr r 


P 


u 


(T — 1 )p 2 




( 6 . 1 ) 

(6.2) 
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where Hq is a predescribed total enthalpy, By sub¬ 
stituting (6.2) into (6.1), a differential equation for p 
is obtained and can be easily solved by numerical in¬ 
tegration. Consequently, p and the entire stationary 
flow field is determined. The solution of this station¬ 
ary vortex can then be superimposed to any uniform 
mean flow with a given streamwise velocity uq. In 
the problem we consider, ug and u r are converted 
to: 

u = 6ye- ar \ v = Sxe- ar \ 

where <5= 0.3 is a given amplitude factor, while p 
and /?, as shown above, are functions of r, with 

r 2 = ((*- u 0 t) 2 + y 2 ). 

This stationary vortex is thus superimposed to the 
background mean flow. The vortex is large and 
strong, since the pressure in the vortex center dips 
down to about 7% of its circumferential value. Due 
to the presence of shocks, the CE/SE scheme is re¬ 
quired to have enough numerical damping. For this 
purpose, a weighted average index a = 2 and nu¬ 
merical dissipation factor e = 0.5 are chosen. The 
boundary conditions are the same as for the ini¬ 
tial oblique shocks computation. We choose a = 
?§ 4 , At = 0.2 and run 900 time steps. Fig.6.10 
demonstrates the interaction process at different 
time t=2,20,38,56,74,92,110,128,146, and 180. The 
vortex propagates downstream while remaining in¬ 
tact before colliding with the oblique shock. With 
further propagation, it begins to interact with the 
first oblique shock ( Fig.6.10 (2-3) ). During the 
interaction, both vortex and shock are deformed or 
distorted. The straight oblique shock first changes to 
S shape and then recovers to its original straightline 
shape. At Stage 2, the collision disrupts the vortex. 
The ruptured vortex can no longer remain intact and 
begins to release its kinetic energy in the form of 
strong (non-linear) acoustic waves (pressure ampli¬ 
tude is as high as 17% of the field maximum). The 
phenomenon is consistent with others’ experiments 
and numerical computations for a normal shock - 
vortex interaction. In Fig. 6.10 (3-4), the disrupted 
eddy continues to emit acoustic waves. Some of the 
waves pass through the reflected shock and are re¬ 
flected from the solid wall. In Fig.6.10 (5), the weak¬ 
ening vortex is flushed further down and collides 
with the second oblique shock. It is then further 
disrupted and releases more energy in the form of 
acoustic waves. At Stages 6,7 and 8, more acous¬ 
tic waves are emitted and the vortex reduces its size 
( to about 1/4 of its original size) and kinetic en¬ 
ergy as well. In Fig.6.10 (9-10), the vortex is flushed 
out from the computational domain along with the 


acoustic waves it generates, and the oblique shocks 
resume their original shapes. 


6.8 Acoustic Pulse / Shock-Wave In- 
teraction 

In order to demonstrate the capability of the CE/SE 
scheme to handle interactions of acoustic waves and 
shock waves, we describe an example of a weak 
acoustic pulse wave passing through a strong shock. 
We use a mesh of 200 x 200 nodes. The domain 
is centered at the origin (0,0), with an extent of 
— 100 < a? < 100 and -100 < y < 100. A steady 
oblique shock at a position along a diagonal of the 
computational domain is precalculated to form part 
of the initial condition for the computation. The 
initial conditions of an isolated acoustic pulse are 
superimposed on this precalculated shock to form 
the initial condition of the given problem. The data 
upstream and downstream of the shock are respec¬ 
tively 

uq = 2.378056, vq = 0, 
po = 1 and po = 0.7142857; 

and 

u 0 = 2.1017481, vq = 0.4062729, 

Po = 1.5807555 and p 0 = 1.3713613. 

A weak acoustic pulse propagating across a strong 
shock is considered. An acoustic pulse, initially cen¬ 
tered at (zo,2/o) = (—75,0), with initial data 

u* = v* = 0, p* = p* = 6e - a [(*-*o) 2 +(y--i/o) 2 ] 

is superimposed on the mean flow, where the ini¬ 
tial pulse amplitude e = 0.001 and a = (ln2)/9. It 
is observed that the oblique shock strength is three 
orders of magnitude larger than the initial ampli¬ 
tude of the acoustic pulse. The pulse propagates in 
all directions with the speed of sound, while being 
carried downstream by the mean flow. During the 
computation, the non-reflecting boundary condition 
described above is enforced at all the four sides of 
the computational domain. 

For such an interaction between a weak ( linear) 
wave and a discontinuous wave, the theoretical exact 
solution is not available. However, the numerical re¬ 
sults obtained with the CE/SE scheme demonstrate 
physically plausible phenomena. Fig. 6.11 illustrates 
the isobars at various time steps. At first, the acous¬ 
tic pulse is blown downstream and propagates freely. 
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As the pulse collides with the strong oblique shock, 
the shock is practically unaffected, while the acous¬ 
tic pulse ring is distorted in its passage through the 
shock, due to different speeds of sound and flow ve¬ 
locities on either side of the oblique shock. 

In other examples described in [9], [15] and [16], 
the interactions of a strong (i.e. nonlinear) acoustic 
pulse, and of weak and strong vortical and entropy 
pulses with a strong shock were computed. Cur¬ 
rently, the CE/SE method is being applied to bench¬ 
mark problems in CAA and to supersonic jet noise 
computations, and has proved to be exceptionally 
accurate. 


6.9 Three-Dimensional Inviscid Flow 
Examples 

The CE/SE 3D Euler solver, which is an extension 
of the ID Euler solver in the same way as is the 2D 
Euler solver, has recently been developed by X.Y. 
Wang and S.C. Chang. See Section 5 for the con¬ 
ceptual description of flux conservation using tetra¬ 
hedrons in 3D. Details of the solver will be presented 
in a paper under preparation [?]. Some preliminary 
numerical results are presented in the following. In 
these results, a structured mesh of tetrahedrons is 
used, with a = 1 and e = 0.5 in the entire domain. 

(a) Oblique shock problem: To help with vali¬ 
dation, the 3D Euler code is used to solve the 2D 
oblique shock reflection problem previously solved 
with the 2D Euler solver in [7]. The structured mesh 
consists of 42x14x14x6 tetrahedral spatial cells in a 
domain 0 < x < 4, 0 < y < 1 , and 0 < z < 1 . 
The steady-state solution obtained with a time step 
At = 0.012 is shown in Fig. 6.12, in which the den¬ 
sity distribution in the entire domain and density 
contours at the plane y = 0 are plotted. 

(b) Flow past a ramp: For further validation, the 
3D Euler code is used to solve a 2D supersonic flow 
°f Moo — 2.2 past a ramp with a compression angle 
0 = 12°, which has previously been solved using the 
2D Euler solver in [22]. 40x10x10x6 tetrahedral cells 
are used in a domain -0.2 < x < 1.8, 0 < y < 

1, and 0 < 2 < 1, and a time step At = 0.005 
is used. The steady-state density distribution and 
density contours at the surface y = 0 are shown in 
Fig. 6.13. 

(c) Implosion/Explosion of a spherical shock wave 
in a cubical box. 41x41x41x6 tetrahedral cells are 


utilized in the computational domain specified by 

^ — — 0<y<2, and 0 < z < 2. The sym¬ 

metry inherent in the problem is utilized, so that 
the initial shock wave configuration is an eighth of a 
sphere, whose geometric center isata; = y = ^ = 0 
and whose radius is 0.8\/3. The interior of the sphere 
is a low pressure region with a pressure ratio of 10 
across the shock. The density and pressure contour 
plots at different time levels obtained using a time 
step At — 0.005 are shown in Fig. 6.14. The shock 
wave, contact surface, and expansion wave in Fig. 
6.14 display physics similar to that seen in the 2- 
D implosion/explosion problem. More detailed de¬ 
scription of the flow physics will be provided in [?]. 


7 Summary and Conclusions 

In the present article, we reviewed the method of 
space-time conservation element and solution ele¬ 
ment (the CE/SE method, for short) for the nu¬ 
merical solution of conservation laws. We described 
several CE/SE schemes for computing fluid flows, 
and touched upon other CE/SE schemes and exten¬ 
sions. Our descriptions emphasized the geometry of 
the space-time discretization. 

An ideal solver for smooth flows must be neutrally 
stable, explicit and two-level, and must be such that 
the discrete equations are invariant under space-time 
inversion. The CE/SE non-dissipative Euler solvers 
for isentropic flows meet all these requirements. In 
the present article, we described the non-dissipative 
ID and 2D Euler solvers in terms of the conserva¬ 
tion of piecewise linear space-time fluxes over dis¬ 
crete space-time volumes. Thus, given the space- 
time discretization, the schemes have a simple spec¬ 
ification in terms of flux conservation. When shock 
waves are present in the solution, numerical dissi¬ 
pation must be introduced into numerical schemes 
in a controllable fashion, to model the irreversibil¬ 
ity in the exact solution. We described the shock¬ 
capturing ID Euler solver, which is a modification 
of the non-dissipative solver. The added numeri¬ 
cal dissipation has a simple geometric description 
and a straightforward generalization to the 2D case. 
The Navier-Stokes solvers, not described here, re¬ 
duce to the non-dissipative solvers when the physi¬ 
cal viscosity vanishes, and hence the latter is never 
overwhelmed by numerical dissipation. 

The key strategies that enable the CE/SE schemes 
to avoid the limitations of the upwind schemes are: 
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(i) The more general form of the conservation laws, 
i.e., the integral form, is cast in a form in which 
space and time are treated on an equal footing. This 
gives flexibility in the shape of the space-time con¬ 
servation elements, which is useful for defining CEs 
when, for e.g., sources are present in the CE. (ii) A 
staggered space-time mesh is employed. This results 
in the simplest stencil. It also obviates the need for 
interpolation of fluxes at the interface between CEs. 
Thus, there is no need for an approximate Riemann 
solver. Hence, characteristics-based upwind-biasing 
methods, which are complicated and strictly valid 
only for smooth solutions, are avoided. There is thus 
also no compromise in the symmetry of treatment 
of the spatial fluxes. This also has implications for 
flows in multiple spatial dimensions. For the com¬ 
putation of such flows, upwind techniques must use 
directional splitting with its attendant difficulties. 
The CE/SE method in multiple spatial dimensions, 
on the other hand, does not involve any directional 
splitting, (iii) The flow property gradient is treated 
as an additional unknown in the CE/SE schemes. 
Therefore, there is no need for reconstruction of the 
flow gradient by polynomial curve fitting over neigh¬ 
boring mesh points, and for the subsequent use of 
complicated flux limiters, (iv) Space-time fluxes are 
conserved at both the local and global level. The 
condition of flux conservation, rather than any ex¬ 
trapolation, links the solution at a mesh point with 
its neighbors at the previous time level. This empha¬ 
sis on the integral conservation law is critical for ac¬ 
curate flow simulations, particularly if they involve 
long marching times and/or regions of rapid change 
(e.g., boundary layers and shocks). 

We reproduced here several numerical results ob¬ 
tained with various CE/SE flow solvers. The re¬ 
sults included a demonstration of extremely simple 
yet highly effective non-reflecting boundary condi¬ 
tions for the extended Sod’s shock-tube problem. 
The CE/SE solver for the scalar convection-diffusion 
equation was shown to be accurate in all Reynolds 
number regimes. The CE/SE solver for the ID and 
2D Euler equations with source terms simulating 
detonation was also shown. We reproduced numer¬ 
ical solutions obtained with the 2D CE/SE Euler 
solver, including the process of diffraction of a shock 
wave around a wedge and the implosion/explosion of 
a polygonal shock wave in a box, as well as compu¬ 
tational aeroacoustic phenomena involving the inter¬ 
action of strong shocks and weak acoustics as well 
as strong shocks and strong vortices. The results 
reproduced here are only some of the difficult prob¬ 
lems readily solved with CE/SE schemes; see [1-31] 


for more examples. 

We remark here that the CE/SE schemes devel¬ 
oped thus far are characterized by simplicity, gen¬ 
erality of applicability and second-order accuracy 
in space and time. The simplest possible stencils 
are employed. The 2D spatial mesh is constructed 
from triangles, and the 3D spatial mesh will be con¬ 
structed from tetrahedra. Triangles and tetrahedra 
are the simplest polytopes in 2D and 3D, respec¬ 
tively. The ID and 2D Euler solvers bear a remark¬ 
able resemblance to the solvers of the ID and 2D 
scalar convection-diffusion equations, respectively, 
with the discrete equations in the former two be¬ 
ing matrix versions of the scalar equations in the 
latter two. All of the above schemes are character¬ 
ized by virtually the same properties. Furthermore, 
the viscous flow solvers are designed to reduce to the 
respective non-dissipative solvers when the physical 
viscosity vanishes. The CE/SE method thus rep¬ 
resents a new unified framework for the numerical 
solution of conservation laws. 
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Figure 2.1 —Space-time conservation elements for methods using a fixed 

spatial domain: (a) one space dimension (b) two space dimensions 



Figure 2.2—A space-time conservation element with an 
arbitrary space-time domain 
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Figure 3.1 —The space-time mesh 
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Figure 3.3 —The CEs associated with the mesh point (j,n) 
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Figure 4.1 — A spatial domain formed from congruent triangles, 
showing the spatial projections of the mesh points 
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Figure 4.3 — The SE associated with the point G' 



Figure 5 .1 — Spatial projection of part of a 3D space-time niesii, 
showing the construction of a CE 
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Fig. 6.6 (a) Front Pressure of Unstable ID ZXD Detonation Wave 




Fig. 6.6 (b) — Pressure Contours of Unstable 2D ZND Detonation Wave 
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Fig. 6.7 Diffraction of a shock wave around a wedge: 
Density contours of the CE/SE solution 
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Fig. 6.8 Density contours for implosion/explosion of a polygonal shock in a square box 
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Fig. 6.9 Density contours for implosion/explosion of a hexagonal shock in a square 
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Fig. 6.11 Isobars for interaction of an acoustic pulse with a shock wave 
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Figure 6.12 The oblique shock problem. 



Figure 6.13 Flow past a 3-D ramp. 





















































(a) t=0.4 (b) t=0.55 




(c) t=0.85 (d)t=1.15 

(a) Density contours 



(c) t=0.85 (d) t=1.15 

(b) Pressure contours 

Figure 6.14 Implosion/explosion of a spherical shock wave at different time levels. 
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